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Abstract We present a double site-bond percolation model to account, on the one hand, for the vascu- 
larization and/or resorption of biomaterial implant in bones and, on the other hand, for its mechanical con- 
tinuity. The transformation of the implant into osseous material, and the dynamical formation/destruction 
of this osseous material is accounted for by creation and destruction of links and sites in two, entangled, 
networks. We identify the relevant parameters to describe the implant and its evolution, and separate 
their biological or chemical origin from their physical one. We classify the various phenomena in the 
two regimes, percolating or non-percolating, of the networks. We present first numerical results in two 
dimensions. 

1 Introduction 

The interest in using natural coral as biomaterial in order to replace missing bone tissue in both or- 
thopaedic and maxillo-cranio-facial surgery has been recognized for many years. This natural material is 
now widely replaced by synthetic one, like hydroxyapathite or bioglasses with various mineral composi- 
tion and porosity. All these materials should replace bony tissue, or favor the formation of bony tissue in 
contact with prothesis, without inapropriate response from the human body (biocompatibility) [1,2,3 . 

In all these processes involving natural or synthetic biomaterials, the ability of biomaterial implants 
to transform into neoformed bones is of particular interest. This ability depends on various parameters 
of the implant. Some of them are of biological or chemical nature, the others are of physical origin. We 
propose in this study a general framework based on percolation theory to clearly identify the relevant 
parameters to describe the structure and the evolution of such biomaterial implants, and the scale, in 
space, they correspond to. 

Our study is in the continuation of both experimental and theoretical previous investigations. On 
the experimental side, the full characterization of the structure and evolution of in vivo coral implants 
have been done using nuclear physics methods - Both neutron activation analyzis and particle 

induced X-ray emission methods allow to determine the elementary composition of implant as a function 
of time. These methods are very sensitive and enable to quantify the trace elements to highlight the 
vascularization process [7| and to observe the bone/biomaterial implant interface. 

In the case of coral implants, these studies have shown the following, very general, pattern: 
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1. during the first few weeks, the mineral composition of the biomaterial implant does not change sig- 
nificantly. This phase corresponds to the invasion of blood cells and to the vascularization of the 
implant; 

2. once this vascularization is achieved, biomaterial resorption and bone remodelling start to operate 
simultaneously. 

On the theoretical side, a first attempt to understand this evolution has been proposed in [5] in the 
framework of the percolation theory. In a simplified model, the implant was characterized by a network 
of sites which can be either "occupied" (or activated) or "unoccupied" (unactivated). The activated sites 
in the initial implant were identified with empty pores, while the unactivated sites were identified with 
biomaterial implant. The initial implant is thus characterized by a single parameter, the probability po 
to have an activated site. A step-by-step process was then proposed to describe first the vascularization 
of the empty pores (which were assumed to be all interconnected) untill the connected vascularized pores 
form a "percolation cluster", i.e. a cluster occupying more than 50 % of the volume of the implant. 

When this happens, we say that the "percolation threshold" has been reached. In the percolation 
cluster, each pore is vascularized and linked to at least one of its nearest neighbor. Once the percolation 
threshold is reached, the osteoblastic affixing responsible for osseous neoformation can start, together with 
remodelling of neoformed bone. This is done in the numerical simulation by giving the possibility, to the 
sites, to transform into a third state associated with the formation of osseous material. The time evolution 
of this transformation is described by a set of parameters, qi, giving the probabilities to transform a site 
from one state (biomaterial, empty or vascularized pores, or osseous material) to another one. 

From a general point of view, the time evolution of implants of biomaterial does depend on its 
physical properties like its porosity, the size of individual pores or the interconnection between them, as 
well as on the chemical and biological properties of the biomaterial itself, like its speed of resorption or its 
bioactivity. In order to account for all these properties, we need a more flexible and general simulation. 
We shall thus extend this simple model in three main directions. 

1. First of all, we shall characterize the physical structure of the implant by giving the possibility to 
each activated site to be connected, or not, by a link to its nearest neighbor. The link between two 
activated sites will be associated to the possibility for the pores to be interconnected. The structure of 
the initial implant is thus characterized by two parameters: the probability, pg, to have an activated 
site, and the probability, q^, to have a link between two neighboring activated sites (empty pores); qo 
will be called the connectivity rate. 

2. The time evolution of the implant will explicitly incorporate the phase space volume available for 
the transformation of sites and links. This is necessary in order to recover the scale invariance of our 
result with respect to the change of the longitudinal dimension of the implant. 

3. Finally, we shall consider an original double percolation lattice in order to be able to account for 
both the vascularization of the implant, forming a blood network, and the formation of a second, 
interconnected, network associated with the osseous material, called osseous network. We can thus 
address the question of the mechanical solidity of the implant, and its evolution in time, as well as 
the compatibility of the two percolation clusters in the blood and osseous networks. 

In a site-bond percolation model, the system is characterized by a network of sites and links 
connecting two neighboring sites, represented on a lattice in D dimensions. The size of the lattice is 
determined by the number, iVj, of sites regularly distributed on each dimension, i, of the lattice. The 
sites can be activated with a probability p, while each pair of neighboring activated sites can be linked 
with a probability q. A cluster is defined by any assembly of activated sites linked together. Percolation 
theory ( f illT(TlllTUT21 tells us that the state of the system can be characterized, in a universal way, by a 
critical line in the two-dimensional (2D) (p, q) diagram above which there is a unique large cluster which 
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occupies most of the size of the system. This cluster is called the percolation cluster. In an infinite - or in 
practice very large - lattice, the transition between the state with (infinitely) many small clusters and the 
one with a single very large cluster is very sharp, which defines unambiguously the critical line and the 
critical points at q = 1 (p®) and at p = 1 (q^)- A typical phase diagram corresponding to this transition 
is shown on figure [l] 
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Fig. 1 Phase diagram of a typical site-bond percolation model in a two dimensional lattice corresponding to a 
disk of diameter L = 800, i.e. which has L sites on its diameter. The region in grey refers to the domain where a 
percolation cluster exists. 

The interest to model the evolution of biomaterial implants within percolation theory is two-fold. 
Firstly, this modeling can be used to identify the relevant parameters to describe the evolution of the 
system, and to clearly separate what is due to the physical properties of the system from what is due 
to its biological or chemical structure. Moreover, it also shows that the main properties of the system 
do depend on the cluster properties of the networks, i.e. whether the system is in the percolating or 
non-percolating phase. Secondly, the modeling of the system is important in order to investigate rather 
easily the influence of physical, biological or chemical properties of the biomaterial implant on its time 
evolution and on the final structure it can lead to when it is embedded into bone tissues. This may avoid 
in-vivo experiments. The consideration of both networks, blood and osseous, is particularly suited to 
address two main properties of neoformed bone: its vascularization and its mechanical strength. 

The plan of the article is the following. We detail in section[2]the construction of the two, blood and 
osseous, networks and present their general static properties. In section [3j we present the time evolution 
of these networks and characterize the dynamical connection between them. We discuss our numerical 
results in a two dimensional simulation in section [4j Conclusions are presented in section [5j 

2 General properties of the double percolation model 

2. 1 Definition of the blood and osseous networks 

The construction of the two lattices is motivated by the structure of the bone tissue in which the bio- 
material is implanted. With a porosity close to 50%, the bone tissue (cortical bone) is composed of two 
entangled networks of similar structure. We thus consider a first lattice of size N. It is assumed to be 
regular for simplicity, with a lattice spacing a. The second lattice is deduced from the first one by a 
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translation in the two dimensional plane by half the lattice spacing in both dimensions. It is of size N — 1 
with the same lattice spacing. Note that in two dimensions, there is no alternative ways to construct 
these two lattices. This is not anymore the case in three dimensions. On the first lattice, we shall define 
the blood network, while the osseous network will be considered in the second, dual, lattice. 

We shall consider in our study a simple 2D geometry for the biomaterial implant. It is assumed to 
be a disk with L sites on the diameter, which is embedded in a piece of bone represented by a square of N 
sites, as indicated in figure [2] for the case of a very small piece of lattice. This geometry is chosen in order 
to mimic the physical case of a cylindrical implant in the transverse plane, as studied experimentally in 
[!J[5l[6] . Various typical configurations of the sites and links are illustrated on figure [2] 
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Fig. 2 left plot: general geometry of the blood and osseous networks, with the initial biomaterial implant in grey, 
labeled by (B), embedded into an osseous network (O); right plot: various states of the sites and links in the two 
networks (see text for more explanations) . The blood network is shown by solid lines for the links and squares for 
the sites, while the osseous network is shown by dashed lines for the links and disks for the sites 



The sites of these networks can be in the following states: 

— For the blood network, the unactivated sites, represented by small full squares, are associated with 
biomaterial implant. The activated sites, represented by large squares, are associated with the absence 
of matter (pores). They may correspond to two different states: either empty pores (large empty 
squares) or vascularized pores (large full squares). 

— For the osseous network, we assume for simplicity that all sites are activated and correspond to matter. 
This matter can be either biomaterial implant (large empty disks) or osseous matter (large full disks). 
We may easily generalize the simulation by allowing the sites in the osseous network to be unactivated. 
This may corrrespond for instance to the presence of defaults in the initial biomaterial implant, or to 
trace elements. 

The links between two neighboring, activated, sites of these networks can be in the following states: 

— For the blood network, the unactivated links, represented by small solid lines in figure [2] are associated 
with matter, either biomaterial implant or osseous material, depending on the nature of the perpendic- 
ular link in the osseous network. The activated links are associated with the absence of matter. They 
may correspond to two different states: either empty, i.e. they correspond to interconnected empty 
pores (double lines) or vascularized (large solid lines), depending on the evolution of the system as 
explained in section [3] 

— For the osseous network, the unactivated links, represented by small dashed lines in figure[2]are associ- 
ated with the absence of matter, i.e. they correspond to either empty or vascularized interconnections 
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depending on the nature of the perpendicular link in the blood network. The activated links are asso- 
ciated with matter. They may correspond to two different states: either biomaterial implant (double 
dashed lines) or osseous material (large dashed lines), depending on the evolution of the system as 
explained in section [3] 

Note that in our 2D lattice, any activated link between two neighboring sites in the osseous network 
should correspond to an unactivated link between to perpendicular neighboring sites in the blood network, 
and vice versa, since these two links cross each other. This will not be the always the case anymore in 
three dimensions. The precise state of the sites and links in both networks will be determined by the 
initial conditions and the evolution process which will be discussed in section [3] 



2.2 The phase diagram in two dimensions 

Before entering into the details of the time evolution of the system, we shall hrst determine the static 
properties of the networks we are considering. We can refer here only to one of the two networks since the 
geometry of both networks is nearly identical. We shall consider the disk geometry for the lattice (inner 
disk labelled B in figure [2]). 

We consider two typical values of the size L of the lattice, L = 40 and L = 800, and three situations: 

(a) The situation where p = 1 and we determine the critical parameter q c = q®. It corresponds to the 
particular case of bond percolation. 

(b) The situation where q = 1 and we determine the critical parameter p c = p®. It corresponds to the 
particular case of site percolation. 

(c) The situation where p — q (first diagonal in the phase diagram) and we determine the critical value 
Pc = q c - 

In order to calculate these critical parameters, we have to follow the size of the largest cluster as a function 
of p and q. For each value of these parameters, we calculate the number of sites, denoted by N s , in the 
largest cluster. The total number of activated sites is, by definition of p, given by pNr, where Nt is the 
total number of sites in the lattice. We can thus calculate the size of the largest cluster with respect to 
the size of the largest possible cluster. It is given by the coefficient c 

C= pNr- (1) 

We repeat the simulation a large number of times and calculate the mean value of c. This value is plotted 
on figure [3] for the three situations mentioned above. Note that we may also consider the number of 
activated links in the largest cluster, or the sum of both sites and links, to define the coefficient c. The 
results are almost identical for all these definitions. 

We see on these figures that the transition between the percolating and non-percolating phases 
is extremely rapid for L = 800, as expected. For small values of L, like L = 40, the transition is still 
abrupt. We shall therefore define the critical parameters p c and q c by the criterium that c should be 
larger than a given value cq- This parameter cq is called the percolation criterium. We shall choose in our 
study Co = 0.4. The results of figure [3] ensure that our numerical simulations will be rather insensitive 
to the choice of the percolation criterium, provided c has no extreme values. The values of the critical 
parameters p®, q® and (p c = q c ) for the three situations (a),(b),(c) are given in table [l] for two typical 
sizes of the lattice. The phase diagram for L — 800 has already been given in figure [T] 
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Fig. 3 Transition between the percolating and non-percolating phase, for L = 40 (dashed line) and L = 800 (solid 
line), for the three situations (a),(b) and (c) discussed in the text. 
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Table 1 Critical parameters p° , g° and (p c , g c ) for the two different lattice sizes, for cq — 0.4. 



2. 3 Application: solidity of the percolation cluster 



The size of the percolation cluster, in any of the two networks, is one of the property of the cluster. 
It however does not give any information on the topology of the cluster. There are several ways to 
characterize this topology p~0llTTTIT2] . Keeping in mind the physical interest of our simulations, we shall 
consider a way which can be interpreted as the solidity of the cluster. 

We shall define the solidity coefficient (number between and 1) as the capability of the percolation 
cluster to persist if we remove a given number of its links. In a given percolation cluster, for a given value 
of (p,q), with Umax activated links, we successively and randomly remove n links. Doing this N s times, 
we call M n the number of times a percolation cluster still persists after the removal of these n links. We 
can thus define the two parameters 



S n (p,q) = ^ }L and 



(2) 



The parameter S n is shown in figure refsolif as a function of q s , for three typical points in the percolating 
phase of the phase diagram and for a lattice size L — 40. The solidity of the cluster can be defined by 
the value of q s , denoted by q°, for which S n = 1/2. These three points are all on the diagonal p = q. The 
first one is just near the transition line. It shows that the percolation cluster is in that case very easily 
broken, as expected, and g° is very close to 0. For p = q = 0.85, we have q® = 0.24. For p and q close to 
1, the cluster is expected to be very stable, and we find q$ — 0.49. Not surprisingly, this value is identical 
to the critical value q® for this size of the lattice, as given in table [l] 
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Fig. 4 Characterization of the solidity of the percolation cluster by the coefficient S n in Q), for a lattice size 
N = 800, and for three typical points (p, q) in the phase diagram indicated directly on the plots. 

3 Double diffusion in the blood and osseous networks 



3.1 Initial conditions 

As already shown in figure[2j the system we consider has two different parts. The first one is the biomaterial 
implant under study (the inner disk, B, on this figure). The second one is the osseous material itself 
surrounding the implant, and denoted by O. 

We have thus to define two types of limiting conditions. The first one is at the external border of 
the osseous material, i.e. the external border of our square lattice. We shall assume in this case that all 
the external links of the two, blood and osseous, lattices are all connected to a (virtual) very large cluster, 
and are thus all activated. The links to this external blood network are assumed to be vascularized, as 
well as the first sites of the blood network they are connected to. The same is true for the osseous network: 
all the links to the external osseous network are assumed to be osseous material, as well as the first sites 
they are connected to. An example of such boundary condition is shown in figure [5] 



Fig. 5 Boundary condition at the external border of the lattice, for the blood network. A similar boundary condition 
is also taken for the osseous network. 



In this way, we construct two initial clusters: the vascularized blood cluster and the osseous cluster. 
Both are supposed to be connected to the entire body in which the biomaterial is implanted through the 
boundary conditions. The evolution of the system will thus start from these two clusters, according to 
the rules we shall detail in the next subsection. 

The second limiting condition we have to deal with is the interface between the biomaterial implant 
and the bone in which it is embedded, as shown by the circle in figure [2j This interface can be modelized 
according to the experimental conditions. In our study, we shall consider a perfect interface, in the sense 
that there will be no discontinuity for the sites: each site will belong either to the osseous part or to the 
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biomaterial part of the system. The links in the blood network will be activated with the probability qq, 
and the non-activated links will correspond, in the osseous network, to biomaterial implant. 

3.2 Evolution of the system 

The evolution of the system is done according to very simple rules. Starting from a randomly chosen site 
in one of the two initial clusters (blood or osseous), we perform a one-step evolution of the cluster by 
activating one of the neighboring sites, and/or its link to the original chosen site, which is not already 
connected to it by a vascularized (ossified) link in the blood (osseous) network. This activation is done 
with a probability qi. This probability depends a priori on the nature of the two neighboring sites, as well 
as on the nature of the link between them. We shall discuss the determination of qi in the next subsection. 
The blood or osseous initial clusters are thus enlarged to include this new site and its link. 

The choice of the site from which the displacement takes place is done in the following way. We 
first calculate the phase space available for the displacement in both networks. We denote this phase 
space by Rb and R for the blood and osseous networks, respectively. The total phase space is thus: 

R = R b + R . (3) 

It is given by the total number of free sites available for a displacement. To calculate Rb and R , we 
classify each site of the two clusters according to the number of free sites it is surrounded by, i.e. the 
number of sites which is not already connected to it by a vascularized (osseous) link in the blood (osseous) 
network. This number can be k = 0, 1, 2 or 3 for a 2D lattice. We show in figure [6] some simple examples 
for k = 2. Each site belongs thus to one among four classes, for the two networks. These classes are 




Fig. 6 Some examples of configuration corresponding to k = 2 for the site shown in the middle, in the blood 
network. 

denoted by C^ . If we call Njf the total number of elements in each class C^ Q , the phase space volume 
i? fc:0 is given by 

3 

R b>0 = Y,kNl . (4) 
fc=i 

In order to restrict ourself to a choice of a random number between and 1, independently of the size 

0kN k 
with Pb — • The 

simulation proceeds through the following successive steps, once the eight classes C\ a have been formed. 

1. We first choose a random number po, between and 1 to determine a class C\ , according to the 
representation of figure [7J The example shown on this figure determines the class C^. 

2. In the chosen class, Cf , we randomly determine a given site i. 

3. From this site, we randomly determine, if necessary, a free neighboring site /. 
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Fig. 7 Representation of the phase space volume in the interval [0,1]. The size of each segment is completely 
arbitrary in this example. 

4. This new site / is thus added with a probability qi to the cluster, and its link to the site i is activated. It 
is "vascularized" if i belongs to the blood network, and "ossified" if it belongs to the osseous network. 

5. The time r is increased by an amount At given by 



— )■ t -\- At , with At = - 



lnjpi] 
R 



(5) 



where pi is a random number between and 1. 
6. The classes C b fe Q are recalculated. 



3. 3 Definition of the parameters 

The dynamics of the evolution of the blood and osseous clusters is entirely determined by the value of 
the parameters qi. These parameters depend on the nature of the connected sites denoted by i and /, 
and on the nature of the link between the two sites. We indicate in table [2] all possible values. 

The parameters q V) qb,q r and q Q stand for the vascularization of empty pores or interconnection 
between empty pores, the bioactivity of the biomaterial, the resorption of biomaterial and the cre- 
ation/destruction of osseous material, respectively. The parameter q v is of physical origin. It depends 
mainly on the size of the pores. In the case of coral implants, these ones are large enough (macrop- 
orosity) to enable the flow of osteoblast/osteoclast cells. The parameters qb and q r refer to properties of 
the biomaterial. They are of biological and chemical nature (chemical properties of the biomaterial and 
biological properties of the blood fluid). The first one, called bioactivity, corresponds to the ability 
of the biomaterial to interact with the biological fluid, and transform into osseous material. The second 
one, q r , called resorption, corresponds to the ability of the biomaterial to be resorbed by the biological 
fluid. The parameter q a refers to the osseous material itself. It depends mainly on the biological fluid and 
corresponds to the ability to create/destroy osseous material. 

The parameters a and (3 are of physical origin. They refer to the amount of matter which should 
be transformed. When both a link and a site has to be transformed, as compared to the transformation of 
a site only, the probabilities qi should be reduced by a factor a, with a < 1. This is the case for instance 
when both an empty pore and an empty interconnection between the two pores should be vascularized 
as compared to the vascularization of an empty interconnection only. The second one is done with a 
probability q v while the first one is done with a probability aq v . It is also the case when both sites and 
link associated with the biomaterial are resorbed and vascularized - with a probability aq r - as compared 
to the transformation of only one link - with a probability q r . 
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The parameter (3 is introduced when the transformation of a link can be achieved from both sites 
connected to it. This happens for instance in the diagram on the left of figure We have therefore 
(3 > 1, with the constraint that the final probability should be less than 1. We have typically j3 = 1.1. 



Blood network Osseous network 



Site i Link Site / qi 


Site i Link Site / qi 


vas. vase. vas. 
int. pore q v 
biom. j3 q r 
bone j3 q 


bone vase. bone /3 q 
int. pore 
biom. 
bone 


vas. vase. pore 

int. pore a q v 
biom. q r 
bone q a 


bone vase. biom. qi, q a 
int. pore 
biom. 
bone 


vas vase. biom. 
int. pore 
biom. a q r 
bone q r q a 





Table 2 Probabilities to diffuse from one site i to another site f , for both the blood and the osseous networks. The 
abbreviation vase., int. pore, and biom. stand for vascularized, interconnected pore and biomaterial respectively. 
See text for the interpretation of a, ft and the parameters qi. When a situation is impossible, it is indicated by —, 
while the parameter qi is zero when there is no transformation possible. 



Since the parameters qi are of biological and/or chemical origin, they may change depending on 
the properties of the biological fluid and on the nature of the biomaterial implant, following the analysis 
of [5]. Before the complete vascularization of the implant, i.e. before the percolation threshold is reached 
in the blood network, the dominant mechanism is the vascularization of the biomaterial implant through 
interconnected empty pores, if the size of the pores is large enough. This is the case for coral implant 
where the mean size of a pore is 150 fim. The dynamics is thus governed by q v , and q r _ are very small. 

Once the percolation threshold is reached in the biomaterial implant, the circulation of blood in 
the blood cluster can provide the necessary cells (osteoblasts and osteoclasts) to resorb and ossify the 
implant. The probabilities q r and q become then sizeable. In order to refer to these two regimes - before 
and after the percolation threshold - we consider both values qf a and q^ al respectively. 

3.4 Scales in time and space 

By construction, our modeling has no intrinsic scale for the physical size of the lattice, as well as for 
the time evolution of the system. The lattice is entirely fixed by the number of sites L, while the time 
r is a dimensionless parameter which evolution is given by ([5]). The physical dimension and time should 
therefore be fixed by identification with some given experimental data. 

As far as the time scale is concerned, we have two typical scales in the evolution process of coral 
implants in bones, as shown in figure [8j The first one is called the percolation time t c . It is the time at 
which the concentration of coral implant starts to decrease significantly. The second one is the half-life 
time ti/2 which characterizes the time at which half of the coral implant has been transformed. For 
convenience, it is measured relatively to t c , i.e. the absolute half-life time is t c + 1 1 / 2 - 
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The general time scale to of our simulation, defined by 

t = t T , (6) 

can be fixed by identifying t c with the time r c at which the percolation cluster has been reached in the 
blood network: 

to--- (7) 

This scale is universal for a given geometry of the lattice. Once this is done, we may thus calculate the 
time fi/2, or in other words, the ratio ti/ 2 /t c which is a dimensionless quantity. This scale being fixed, 
we can make predictions for any size, geometry and interface of the implant of a given biomaterial. 

The scale in space is fixed from the geometrical structure of the biomaterial implant. The lattice 
spacing a should be identified with the size of the relevant structures in the biomaterial implant. This is 
neither the overall size of the implant, given by La, nor the microscopic size at the level of the molecular 
structure of the implant. It is given by the size of the typical structures of the biomaterial at the mesoscopic 
scale, like for instance the average distance between the pores of the implant. In the example we choose 
in this study, it is of the order of the typical size of the pores of the coral implant, of about 150 \xm. 
Since the implant has a transverse size of about 6 mm in the experimental conditions of [6] , we choose a 
lattice size of L = 30. 

4 Numerical results 

4--1 Choice of the parameters 

The values of the parameters we choose in our simulation are given in table [3] These parameters refer 



qv qb 


qr qf 


qr Qo 


a 


P 


1 0.2 





0.4 0.7 


0.5 


1.1 



Table 3 Values of the parameters we use in our simulation. 



qualitatively to the transformation of coral implant, with pores of large size so that the vascularization 
process dominates. We therefore can take as a reference q v = 1 in order to speed up the numerical 
simulation. The effect of resorption of biomaterial implant before the percolation threshold has been 
reached in the blood network is taken to be negligible compared to the vascularization, and we take 
qf- Q = 0. The values of q> govern the time evolution of the biomaterial implant once the percolation 
threshold has been reached in the blood network, as shown in figure [8| 

The physical parameter a is taken at a typical value. If the pores and the interconnection between 
the pores are of similar mean sizes, we can take a — 0.5. The parameter fj is fixed to 1.1. 

4-2 Results of dynamical simulations in two dimensions 

In order to remove any sensitivity to initial boundary conditions at the interface between the biomaterial 
implant and the bone, we shall proceed in two steps. We first prepare the formation of the bone sur- 
rounding the biomaterial implant (surface denoted by O in figure [2]), by an (accelerated) evolution of the 
system with an arbitrary structure with initial parameters Po,Qo chosen in the percolating phase, and 
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with all qi equal to 1. Once this is achieved, we start the transformation of the biomaterial implant itself, 
i.e. the inner disk (_B) in figure [2j with physical parameters chosen according to our discussion in section 



3.3[ and given in table [3] The percolation criterium will then be defined within the biomaterial part of 
the simulation only. 



Each of these two simulations is done according to the rules detailed in section 3.2 At each time 
step, we know the number of sites and links of each kind. We can thus calculate the densities of remaining 
biomaterial or of neoformed bones within the domain B of figure [2] To do that, we however have to specify 
the physical sizes associated to sites and links. In this study, we consider a simple configuration where the 
mean size of the pores is equal to the mean size of the links. The density of biomaterial, denoted by pbio, 
and the density of neoformed bones, denoted by pbon, and measured relatively to the initial (subscript 0) 
or final (subscript oo) states respectively, are thus given by 



Pbio 



Pbon 



"bio 



%on 



"■bio 



)o ' 



(n 



bon 



(8) 
(9) 



where n s k {n l k ) is the number of sites (links) of nature k, with k = bio or k = bon for biomaterial or 
neoformed bone, respectively. We stop the simulation when the densities do not change by more than 
1%. We repeat our simulation a large number of times (about 100), and take the mean value of the 
calculated densities. These mean values are shown in figure [8] for a typical lattice of L = 30 sites and 
for a percolation criterium of c a = 0.4, with the parameters indicated in table [3j These parameters have 
been chosen to reproduce the experimental results of coral implants in porcine cortical bones [5]. 




2 4 6 8 10 12 
t(weeks) 

Fig. 8 Evolution of densities of biomaterial and neoformed bone in porcine cortical bones, as a function of time. 



4-.S Trajectories in the phase diagram 
The evolution of the system can be characterized by the time evolution of the various densities defined 



i|9 ). It can also be characterized by its trajectory in the phase space diagram of figure [T] The initial 
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composition of the biomaterial implant is specified by a point (p ,Qo) in the phase diagram. At regular 
time intervals, we calculate the probability p that a site in the blood network is an activated site and the 
probability q to have an activated link between two activated sites. We can thus follow the trajectory of 
the system in the phase space diagram. It is indicated in figure [9] for four different initial conditions. 



1.0 

0.8 
0.6 
0.4 
0.2 
0.0 

0.0 0.2 0.4 0.6 0.8 1.0 
P 

Fig. 9 Phase diagram trajectories for four different starting points. 

When we start from initial points in the non-percolating phase, the system evolves until it reaches 
the percolation transition line, and thus follows closely this line to the point (p = 1,9°). Note that for 
this evolution to take place, the value of q< should not be zero, otherwise the system can not evolve. 
When one starts from an initial point in the percolating phase, the system evolves directly to the same 
point, by first increasing the number of activated sites, and then decreasing the number of activated links 
by the ossification process. 

It is remarkable that all trajectories lead to the same fixed point defined by (p = l,q = q®), 
independently of the initial conditions. This is a direct consequence of the continuous dynamical rear- 
rangement of sites and links between the two networks. This, in turn, is reminiscent of the physical 
processes at work in the biomaterial implant and bone tissue with a continuous creation/destruction of 
bone tissue from osteoblasts and osteoclasts cells. 



5 Conclusions and perspectives 

Our modeling, with two entangled networks, is rich and flexible enough to account for many of the 
processes occuring at the mesoscopic scale in the evolution of implants of biomaterial to osseous material. 
These processes can have either a physical, biological or chemical origin. We choose in this study to 
consider the case of coral implants which has a very similar geometrical structure to the (cortical) bone 
itself. 

The two-dimensional structure we consider is of course a simplified situation. It however may 
correspond to a cylindrical geometry in which the longitudinal scale is much larger than the transverse 
scale. It is also much more constraint than a 3D modeling in the sense that the coexistence of the blood 
and osseous networks is very difficult to achieve in the 2D case. In a static simulation, this coexistence is 
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only achieved in a very limited region of phase space near the percolation transition line. In the dynamical 
simulation howether, with a continuous creation/destruction of links and sites in both blood and osseous 
networks, the final configuration is always a coexistence of both networks. This feature is in agreement 
with the dynamical process at work in the bone tissues. 

The extension of our simulation in three dimensions can be done rather easily. The relative position 
of the blood and osseous lattices is however more flexible then in two dimensions. The translation between 
these two lattices can be done in fact in three different ways. This may lead to systems of various com- 
positions: completely homogeneous, with a one dimensional, linear, structure, or with a two dimensional, 
layer, structure. 
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